Cooling of a lattice granular fluid as an ordering process 
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We present a new microscopic model of granular medium to study tlie role of dynamical correlations 
and the onset of spatial order induced by the inelasticity of the interactions. In spite of its simplicity, 
it features several different aspects of the rich phenomenology observed in granular materials and 
allows to make contact with other topics of statistical mechanics such as diffusion processes, domain 
growth, persistence, aging phenomena. Interestingly, while local observables being controlled by 
the largest wavelength fluctuations seem to suggest a purely diffusive behavior, the formation of 
spatially extended structures and topological defects, such as vortices and shocks, reveals a more 
complex scenario. 
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The effort devoted during the last decades to investi- 
gate off equilibrium systems has achieved a series of suc- 
cesses by virtue of a combination of experiments, numer- 
ical simulations, exact theoretical results and clever phe- 
nomenological arguments, but our comprehension of the 
area is far from complete. Among these systems granular 
materials (GM), i.e. assemblies of macroscopic particles 
which dissipate their energy through inelastic collisions 
and frictional forces, have acquired a special rank due to 
their complex phenomenology often intriguing and only 
partly understood jl]. 

The existing theoretical approaches are based either on 
realistic descriptions of the grains or on idealized mod- 
els, which in virtue of their major simplicity lend them- 
selves to analytic solutions or to efficient computations. 
Furthermore such an idealized modeling can deliberately 
sacrifice likeness to reality, in order to identify the phys- 
ical mechanism responsible for the salient features of the 
system. In such a spirit, we shall discuss a minimal model 
based on the simplest rule which describes the inelasticity 
of collisions and local momentum conservation. 

Years ago, S. Ulam showed that an ensemble of elastic 
particles, starting from an arbitrary configuration, con- 
verged to a Maxwell equilibrium distribution, postulat- 
ing a simple redistribution law of the kinetic energies of 
randomly selected pairs to simulate the effect of binary 
collisions in an elastic gas Ben Naim and Kaprivski 
(BK), recently, performed a variation over this theme, 
by letting the particles endowed with a scalar velocity 
to dissipate inelastically a fraction of the relative kinetic 
energy at each collision and found that the total kinetic 
energy decreases exponentially with time |^] . Both mod- 
els fulfill Boltzmann's molecular chaos hypothesis, and 
consequently rule out the formation of dynamical corre- 
lations. On the other hand computer simulations have 
shown the appearance of a shear instability, i.e. vortices, 
during the cooling process, before the spontaneous forma- 
tion of density clusters. In an interesting series of papers, 
Ernst and collaborators have put forward a mesoscopic 



theory of these phenomena, making a connection with 
phase ordering kinetics 1^ . 

In the present work we shall introduce and study a new 
microscopic model which preserves the simplicity of the 
approach d la Ulam, and displays a complexity similar to 
that observed in granular systems. The focus of our study 
will be on the statistics of the velocity field and on its 
spatial and temporal correlations, stressing the analogies 
and the differences with related models aimed to describe 
off equilibrium systems. 

We introduce our dynamical model by associating a 
d-dimensional velocity field Vi with each node of a d- 
dimensional lattice; at each time step a nearest neighbor 
pair j) is randomly selected and the two velocities are 
updated according to the rule: 

vi = Vj + e(-(vi - vj)a)^^((vi - Vj)a)a 

vj = v; - e(-(vi - vj)a)i±^((vi - Vj)<7)a (1) 

where a represents the unit vector pointing from site i 
to j, G is the Heaviside function which enforces the kine- 
matic constraint and a the normal restitution parameter. 
We shall measure time in the non dimensional number 
of collisions per particle. In each elementary collision 
(see Eq.(^) the total linear and angular momentum are 
conserved, whereas a fraction (1 — a^)/4 of the relative 
kinetic energy is dissipated. The inelasticity of the colli- 
sions has the effect of reducing the quantity |(vi — Vj)(T|, 
i.e. induce a partial allignement of the velocities. Here- 
after we report the study performed in two dimensions 
on a triangular lattice |^ . 

The freely cooling process exhibits striking similari- 
ties with the quench from an initially stable disordered 
phase to a low temperature phase in a magnetic system: 
whereas in a standard quench process pO[ | one considers 
the process by which a thermodynamic system, brought 
out of equilibrium by a sudden change of an external con- 
straint, such as temperature or pressure, finds its new 
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FIG. 1. Energy decay for a = 0.9 and a = 0.2 (1024^ 
sites). The dotted line ~ 1/t is a guide to the eye for the 
asymptotic energy decay. In the insets we reported the scale 
dependent temperature, Ta, defined in the text, as function of 
the coarse graining size a for t = 47 and t = 1720. The total 
energy per particle and To- remain nearly indistinguishable in 
the early incoherent regime, but for a < L(t) the thermal 
energy becomes much smaller than the kinetic energy, a clear 
indication of the onset of macroscopic spatial order. 



equilibrium state, in a GM one wants to study the re- 
laxation of a fluidized state, after the external driving 
force (whose role is to reinject the energy dissipated by 
the collisions, keeping the system in a statistically steady 
state) is switched-off abruptly at some time t — Q. The 
rotational symmetry of the order parameter Vi and the 
momentum conserving interaction determine the pres- 
ence of many configurations having comparable dissipa- 
tion rates. Due to their competition the system does 
not relax immediately towards a motionless state, but 
displays a phenomenology similar to that observed in a 
coarsening process. 

One sees from Fig.^ that during the initial stage, the 
total energy per particle e{t) — Yl^ii^Y / ^ is dissipated 
at an exponential rate = (1 — q:^)/4. This can be 
deduced from Eq.(|l|) imposing that each 'spin' fluctuates 
independently of the others. For times larger than t ^ tc, 
of the order of r, the system displays a crossover to a 
different regime, where the cooperative effects become 
dominant and the average energy per particle decay as 
e(t) ~ t~^. Such a behavior agrees with inelastic hard 
spheres simulations (IHSS) reported in [0. As shown 
below, the crossover from one regime to the other is due 
to the formation of a macroscopic velocity field. This 
is analogous to the formation of magnetic domains in 
standard quench processes. After the formation stage 
these regions start to compete to homogenize, causing a 
conversion of kinetic energy into heat by viscous heating, 
i.e. act against the coUisional cooling and lead to a slower 
decay of the energy Q . 

Within the early regime the velocity distribution devi- 
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FIG. 2. Data collapse of the transverse (S') and longitu- 
dinal {S') structure functions for q = 0. and a — 0.9 (system 
size 1024^ sites, times ranging from t = 500 to t = 10*). The 
wavenumber k has been multiplied by Notice the presence 
of the plateaux for the more elastic system. For comparison 
we have drawn the laws x~'^ and exp(— a;^). 



ates sensibly from a Maxwell distribution (corresponding 
to the same average kinetic energy), but displays fatter 
tails, a phenomenon which mirrors the behavior of the 
BK [H model. The existence of these tails seems to be 
due to the lack of spatial correlations, intrinsically ab- 
sent at all times in their model, whereas negligible in 
ours up to tc- When the energy begins to decay as 
the velocity distribution turns Gaussian. Interestingly, 
the smaller the inelasticity, the faster is the energy dissi- 
pation, a phenomenon observed in IHSS 

The most relevant information about the spatial or- 
dering process is contained in the equal-time structure 
functions, i.e. the Fourier transforms of the velocity cor- 
relation function: 



k 



(k,t)v*''(-k,t) 



where the superscripts t, I indicate the transverse and 
longitudinal components of the field with respect to the 
wave vector k and the sum X^fe O'^^^ ^ circular shell 
of radius k. Such structure factors, if rewritten in terms 
of the variable (fct^/^) display fairly good data collapse, 
apart from the large k region, and identify two grow- 
ing lengths L^'^t) (see Figj2|). Considering the sum rule 
e{t) = ^)] we observed that in the early 

'exponential' regime the contribution from the two terms 
is approximately equal, whereas for times larger than tc 
and a not too small most of the kinetic energy remains 
stored in the transverse field, while the longitudinal com- 
ponent decays faster, with an apparent exponent t~^. 

The findings concerning the energy decay, the distri- 
bution of the velocity field, and the growth of L*''(i), 
lead to the conclusion, that, if the observation time is 
longer than the time between two collisions and if the 
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spatial scale is larger than the lattice spacing, the sys- 
tem behaves as if its evolution were governed by a dif- 
fusive dynamics [^. To be more precise, let us con- 
sider a vector field (j>{x, t) which evolves according to 
the law dt4> = starting from a random uncorre- 

lated initial condition. The explicit solution shows that 
4>{x,t) is asymptotically Gaussian distributed, with a 
variance < 4){x,t)(j){x,t) >(x t~'^/'^. The structure fac- 
tors S^''-{k,t) assume a scaling form S^'\k,t) = s{kL{t)) 
where L{t) — \fi. Furthermore, we compared the two- 
time autocorrelation C(t\^t-2) = X]j '^i(^i)'^i(^2)/^ with 
C^{ti,t2) =< 4'{x,ti)(j){x,t2) >, whose expression reads: 



During a short time transient, the autocorrelation 
function of our model differs from C^, since it depends 
on ti — t2, i.e. it is time translational invariant (TTI). 
Later, C{ti,t2) reaches the "aging" regime and depends 
only on the ratio x — ^1/^2- Something similar occurs 
in a coarsening process, where the autocorrelation of the 
local magnetization a{tw, tw +''■) reaches, for large r (but 
T << itu), a constant value m1^^{T), that is the square 
of the equilibrium magnetization. Obviously, for T — s- 0, 
— > 1 and the TTY transient regime disappears. The 
short time transient in our model is analogous to such 
a TTI regime, with the difference that the cooling pro- 
cess imposes a decreasing temperature T{tw) 0, that 
progressively erodes the TTI regime. The same depen- 
dence on the TTI manifests itself in the angular auto- 
correlation: 



N ^ 



i{t + tw) — 6i{tw))- 



Again, for large waiting times tw this function assume 
the diffusive t/tw scaling form, but for a small fixed tw, 
displays a minimum and a small peak before decreasing 
at larger t (see FigJ|). The non-monotonic behavior of 
A{t,tw) suggests that the initial direction of the veloc- 
ity induces a change in the velocities of the surrounding 
particles, which in turn generates, through a sequence of 
correlated collisions, a kind of retarded field oriented as 
the initial velocity. As tw increases the maximum is less 
and less pronounced. 

In spite of these first results, that seem to give sup- 
port to the idea that the model dynamics is purely diffu- 
sive PI , the model is more complex. The main evidence 
stems from the following facts: 

i) the structure functions do not have the typical Gaus- 
sian tails of a diffusive system, due to the non linearity 
represented by the kinematic factor in Eq.(|l|) and the 
shapes of S*'''{k, t) display three different regions: a long- 
wavelength region which is diffusive in character; an in- 
termediate region where the structure functions decay as 
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FIG. 3. Angular autocorrelation function A{t,tw) for dif- 
ferent values of the waiting time t^j and a = 0.9 (1024^ sites). 
The graph on the left shows the convergence to the t/t^ dif- 
fusive scaling regime, for large t^. For small tm, a local min- 
imum is visible (for such a quasi elastic dynamics). In the 
graph on the right the same data are plotted vs t ~ t^'- note 
that the small curves tend to collapse. For higher the 
position of the local minimum does not move sensibly, but its 
value grows and goes to 1 for large t^. 

k~^ with /3 ~ 4; a plateau region where S"*'' decay in 
time with a power law t~^ but remain nearly constant 
with respect to k (for r > 0); 

ii) the Fourier modes interact and an initial shear state, 
obtained assuming the initial configuration to be a plane 
wave, decays into shorter wavelength modes by a mech- 
anism of period doubling; that is to say, contrary to the 
diffusion, plane waves are not eigenmodes. 

The existence of the quasi-elastic plateaux is the fin- 
gerprint of localized fluctuations which, for small inelas- 
ticity, propagate and are damped less than exponentially. 
A small a determines a rapid locking of the velocities of 
neighboring elements to a common value, while in the 
case of a — > 1, short range small amplitude disorder per- 
sists within the domains, breaking simple scaling of 5*'' 
for large k and having the effect of a self induced noise. 

One can characterize such internal noise by means of 
an average local granular temperature To-, i.e. a mea- 
sure of the variance of Vi with respect to the local av- 
erage of V within a region of linear size a. Obviously, 
since when u ^ 00 the local average tends to the global 
(zero) momentum, then —^ e, as shown in the insets 
of Fig.^ For cr < L{t), instead, T^ < e. For quasi elas- 
tic systems T^ exhibits a plateau for 1 ^ ct <C i(t) that 
identify the strength of the internal noise. The local tem- 
perature ceases to be well defined for smaller a due to 
the absence of scale separation between microscopic and 
macroscopic fluctuations in the strongly inelastic regime 
The existence of a L~^{t)k^^ region in the struc- 
ture functions is consistent with Porod's law |l^] and is 
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Velocity gradient 

FIG. 4. Probability densities of the longitudinal and 
transverse velocity increments. The main figure shows the 
p.d.f. of the velocity gradients (-R = 1). The inset shows the 
Gaussian shape measured for 7? = 40 (larger than L{t) for 
this simulation: a = 0.2, t = 620, system size 2048^). 

the signature of the presence of vortices, a sahent fea- 
ture of the cooling process. Vortices form spontaneously 
and represent the boundaries between regions which se- 
lected different orientations of the velocities during the 
quench and are an unavoidable consequence of the con- 
servation laws which forbid the formation of a single do- 
main. With the random initial conditions adopted, vor- 
tices are born at the smallest scales and subsequently 
grow in size by pair annihilation, conserving the total 
charge. By locating the vortex cores, we measured the 
vortex density Pv{t), which represents an independent 
measure of the domain growth, and in fact it decays 
asymptotically {t 3> tc) as an inverse power of time, i.e. 
Ly{t) = pv^^'^ (X £^1"^. The vortex distribution turns 
out to be not uniform for a not too small. Its inhomo- 
geneity is characterized by the correlation dimension d2 : 
H{R) = e(i? - (ri - rj)) - where the r, 

are the core locations. For a ^ \ the vortices are clus- 
terized < 2) i.e. do not fill homogeneously the space, 
whereas at smaller a their distribution is homogeneous 
id2 ^ 2). 

Vortices are not the only topological defects of the ve- 
locity fields. In fact we observe shocks, similarly to recent 
experiments in rapid granular flows p^ . Shocks have a 
major influence on the statistics of velocity field, i.e. on 
the probability distributions of the velocity increments. 
The probability density function (p.d.f.) of the longitu- 
dinal increment 

A,(R) = l^(v.+K-vO-| 

is shown in Fig.^ for R — I (longitudinal velocity gradi- 
ent) in the main frame, and for i? = 40 > L{t) in the 
inset. For small R <^ L{t) the longitudinal increment 
p.d.f. is skewed with an important positive tail, whereas 



for R ^ L{t) it turns Gaussian. The distribution of 
transverse increments (vi+R — Vi) x R, instead, is always 
symmetric, but non Gaussian distributed for small R. A 
similar situation exists in fully developed turbulence [ p^ . 

To conclude, our model provides a link between the 
microscopic rules of granular dynamics and its hydrody- 
namical description. It allows to follow the cooling of a 
granular material and the build up of velocity correla- 
tions, by means of efficient numerical measures of struc- 
ture factors, two-time correlations and topological de- 
fects. The data analysis reveals the presence of vortices, 
shocks and internal noise and suggests the existence of a 
scale separation only in the case of quasi elastic systems, 
which is instead suppressed for large inelasticities. 

Even independently from the problem of granular 
flows, the model represents a simple but unusual phase 
ordering system. In fact, despite the apparently purely 
diffusive aspects shown by one-point quantities, it dis- 
plays anomalous statistics of spatial properties for the 
order parameter field as witnessed by the velocity gradi- 
ent p.d.f. and by the structure functions. 
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